Trabajo Final - Estadistica Espacial

Instalamos los paquetes necesarios

library(dplyr)
library(sf)

Cargamos los datasets

estaciones <- read.csv("data_smn/preprocessed/estaciones_smn_v2.csv")
horarios <- read.csv("data_smn/preprocessed/datohorario20210420.csv")

Observamos como esta compuesto el dataset

glimpse(estaciones)

Rows: 123
Columns: 9
$ NOMBRE           <chr> "BASE BELGRANO II", "BASE CARLINI (EX JUBANY)", "BASE ESPERANZA", "BASE MARAMBIO", "BASE ORCADAS", "BASE SAN MARTI…
$ PROVINCIA        <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "BUENOS AIRES", "BUENOS AIRES", "BUE…
$ LATITUD_GRADOS   <int> 77, 62, 63, 64, 60, 68, 36, 38, 37, 36, 34, 38, 37, 36, 34, 34, 34, 34, 36, 37, 34, 34, 34, 35, 36, 35, 37, 35, 34…
$ LATITUD_MINUTOS  <int> 52, 14, 24, 14, 45, 8, 50, 44, 43, 12, 32, 1, 26, 21, 36, 49, 33, 58, 2, 56, 33, 41, 40, 27, 53, 52, 36, 22, 27, 1…
$ LONGITUD_GRADOS  <int> 34, 58, 56, 56, 44, 67, 59, 62, 59, 61, 58, 61, 61, 57, 58, 58, 60, 57, 59, 57, 58, 58, 58, 60, 60, 61, 62, 57, 58…
$ LONGITUD_MINUTOS <int> 34, 38, 59, 43, 43, 8, 53, 10, 47, 4, 40, 20, 53, 44, 36, 32, 55, 54, 8, 35, 49, 44, 38, 53, 13, 54, 23, 17, 35, 1…
$ ALTURA           <int> 256, 11, 24, 198, 12, 7, 147, 83, 207, 94, 26, 247, 233, 9, 12, 20, 81, 23, 36, 21, 32, 36, 24, 76, 166, 87, 304, …
$ NRO              <int> 89034, 89053, 88963, 89055, 88968, 89066, 87641, 87750, 87649, 87640, 87570, 87683, 87637, 87648, 87571, 87576, 87…
$ NroOACI          <chr> "SAYB", "SAYJ", "SAYE", "SAWB", "SAYO", "SAYS", "SAZA", "SAZB", "SAZJ", "SAZI", "SADO", "SABP", "SAZC", "SAZD", "S…

glimpse(horarios)

Rows: 2,041
Columns: 8
$ FECHA  <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 2004…
$ HORA   <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, …
$ TEMP   <dbl> 20.7, 19.8, 19.5, 19.2, 19.2, 18.9, 19.0, 18.8, 19.0, 19.1, 19.9, 20.3, 21.8, 22.9, 23.2, 23.7, 23.2, 23.4, 22.5, 22.5, 22.1…
$ HUM    <int> 70, 79, 79, 81, 81, 84, 84, 84, 84, 84, 81, 82, 68, 72, 68, 66, 71, 68, 75, 72, 73, 68, 68, 68, 90, 88, 91, 90, 93, 94, 98, …
$ PNM    <dbl> 1020.8, 1020.9, 1020.9, 1020.4, 1020.0, 1020.3, 1020.7, 1021.1, 1021.6, 1022.2, 1022.4, 1021.9, 1021.7, 1021.0, 1020.1, 1019…
$ DD     <int> 70, 70, 70, 70, 70, 50, 70, 50, 50, 50, 50, 50, 50, 70, 70, 70, 70, 70, 90, 90, 90, 90, 70, 70, 50, 50, 50, 20, 50, 20, 50, …
$ FF     <chr> "28", "28", "26", "26", "22", "22", "17", "17", "15", "19", "20", "24", "20", "19", "20", "17", "15", "15", "19", "20", "22"…
$ NOMBRE <chr> "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQU…

Observamos que la variable que representa la velocidad del viento en km/h es de tipo Char, por lo cual debemos convertila en int.

horarios$FF <- as.numeric(horarios$FF)

NAs introduced by coercion

summary(estaciones$ALTURA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    3.0    53.5   145.0   331.5   455.0  3459.0 

summary(horarios$HORA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    7.00   12.00   12.16   18.00   23.00 

summary(horarios$TEMP)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 -20.00   16.10   20.20   19.09   23.80   33.60       1 

summary(horarios$HUM)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  18.00   59.00   72.00   71.18   85.00  100.00       4 

summary(horarios$PNM)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  960.8  1011.0  1015.1  1075.3  1019.2  3133.0      23 

summary(horarios$DD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0    50.0    70.0   118.9   160.0   990.0       1 

summary(horarios$FF)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   0.00    9.00   15.00   16.06   22.00  107.00       1 

hist(estaciones$ALTURA, main = "Histograma de la altura", xlab = "Altura")

hist(horarios$HORA, main = "Histograma de horarios", xlab = "Horarios")

hist(horarios$TEMP, main = "Histograma de temperatura", xlab = "Temperatura")

hist(horarios$HUM, main = "Histograma de humedad", xlab = "Humedad")

hist(horarios$PNM, main = "Histograma de presion atmosferica", xlab = "Presion Atmosferica")

hist(horarios$DD, main = "Histograma de la direccion del viento", xlab = "Direccion del viento")

hist(horarios$FF, main = "Histograma de la velocidad del viento", xlab = "Velocidad del viento")

Creamos unos boxplot para visualizar la distribucionde datos que potencialmente nos interecen para proceder con el analisis

boxplot(estaciones$ALTURA, main = "Boxplot de la altura", xlab = "Altura")

boxplot(horarios$HORA, main = "Boxplot de horarios", xlab = "Horarios")

boxplot(horarios$TEMP, main = "Boxplot de temperatura", xlab = "Temperatura")

boxplot(horarios$HUM, main = "Boxplot de humedad", xlab = "Humedad")

boxplot(horarios$PNM, main = "Boxplot de presion atmosferica", xlab = "Presion Atmosferica")

boxplot(horarios$DD, main = "Boxplot de la direccion del viento", xlab = "Direccion del viento")

boxplot(horarios$FF, main = "Boxplot de la velocidad del viento", xlab = "Velocidad del viento")

library(DataExplorer)
plot_boxplot(horarios, by = "HORA")

plot_boxplot(horarios, by = "TEMP")

Ahora, veamos que tan correlacionadas estan estas variables

library(ggcorrplot)

corr <- cor(horarios[, c(2,3,4,5,6,7)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
 lab=TRUE,
 ggtheme = ggplot2::theme_gray,
 colors = c("#6D9EC1", "white", "#E46726"))

Como correlacionan estas variables con la altura?

Hay datos nulos en estos datasets?

print("Nulos en horarios")

[1] "Nulos en horarios"

which(is.na(horarios))

 [1]  4853  6710  6894  7999  8000  8401  8935  9137  9138  9139  9140  9141  9142  9143  9144  9556  9557  9558  9559  9560  9561  9562
[23]  9563  9564  9565  9566  9567  9605 10976 13017

print("Nulos en estaciones")

[1] "Nulos en estaciones"

which(is.na(estaciones))

integer(0)

Borramos los nulos del dataset

estaciones = na.omit(estaciones)
horarios = na.omit(horarios)

# Analisis de simetria
library(psych)
skew(horarios$FF)

[1] 1.591713

kurtosi(horarios$FF)

[1] 7.23501

Necesitamos convertir las latitudes y longitudes a

library(biogeo)

latitud <- c()
longitud <- c()

for(i in 1:nrow(estaciones)) {
     latitud[i] <- dms2dd(dd = estaciones[i, "LATITUD_GRADOS"], mm = estaciones[i, "LATITUD_MINUTOS"], ss = 0, ns = "S")
     longitud[i] <- dms2dd(dd = estaciones[i, "LONGITUD_GRADOS"], mm = estaciones[i, "LONGITUD_MINUTOS"], ss = 0, ns = "S")
}

estaciones['LATITUD'] <- latitud
estaciones['LONGITUD'] <- longitud

Unimos las dos tablas usando la variable NOMBRE como punto para combinar los datasets

data <- full_join(estaciones, horarios, by = c("NOMBRE" = "NOMBRE"))

glimpse(data)

Rows: 2,040
Columns: 18
$ NOMBRE           <chr> "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO…
$ PROVINCIA        <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA…
$ LATITUD_GRADOS   <int> 77, 77, 77, 77, 77, 77, 77, 77, 62, 62, 62, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 63, 63, 64, 64, 64, 64, 64…
$ LATITUD_MINUTOS  <int> 52, 52, 52, 52, 52, 52, 52, 52, 14, 14, 14, 14, 14, 14, 14, 14, 24, 24, 24, 24, 24, 24, 24, 24, 14, 14, 14, 14, 14…
$ LONGITUD_GRADOS  <int> 34, 34, 34, 34, 34, 34, 34, 34, 58, 58, 58, 58, 58, 58, 58, 58, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56, 56…
$ LONGITUD_MINUTOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 38, 38, 38, 38, 38, 38, 38, 38, 59, 59, 59, 59, 59, 59, 59, 59, 43, 43, 43, 43, 43…
$ ALTURA           <int> 256, 256, 256, 256, 256, 256, 256, 256, 11, 11, 11, 11, 11, 11, 11, 11, 24, 24, 24, 24, 24, 24, 24, 24, 198, 198, …
$ NRO              <int> 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89053, 89053, 89053, 89053, 89053, 89053, 89053, 89053, 88…
$ NroOACI          <chr> "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "SAYJ", "S…
$ LATITUD          <dbl> -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -62.23333, -62.23333, -62.…
$ LONGITUD         <dbl> -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -58.63333, -58.63333, -58.…
$ FECHA            <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 2004…
$ HORA             <int> 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, …
$ TEMP             <dbl> -20.0, -16.2, -16.9, -13.0, -15.1, -16.6, -15.3, -16.8, -0.1, -0.5, -0.9, -0.5, 0.3, 1.1, 1.3, 2.5, 3.9, 3.1, 3.5,…
$ HUM              <int> 70, 72, 80, 73, 86, 66, 77, 85, 88, 82, 77, 75, 84, 92, 99, 91, 43, 49, 67, 81, 67, 62, 77, 83, 74, 78, 75, 75, 74…
$ PNM              <dbl> 975.7, 970.4, 966.7, 964.2, 961.4, 963.0, 964.8, 967.1, 992.9, 993.0, 992.9, 994.6, 995.3, 995.6, 994.9, 992.3, 98…
$ DD               <int> 140, 140, 110, 140, 160, 160, 90, 110, 230, 230, 270, 290, 290, 320, 320, 320, 160, 200, 50, 200, 110, 250, 250, 2…
$ FF               <dbl> 13, 30, 50, 52, 74, 61, 44, 13, 15, 15, 26, 28, 41, 37, 50, 74, 20, 7, 13, 7, 11, 46, 56, 74, 19, 20, 19, 20, 15, …

Eliminamos las observaciones nulas en el dataset

data = na.omit(data)

corr <- cor(data[, c(7,13,14,15,16,17,18)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
 lab=TRUE,
 ggtheme = ggplot2::theme_gray,
 colors = c("#6D9EC1", "white", "#E46726"))

Agregamos los datos por nombre calculando el promedio y desvio del viento

data_agg = data %>%
  group_by(NOMBRE) %>%
  summarise(MEAN_VIENTO_KMH = mean(FF), SD_VIENTO_KMH = sd(FF), LONGITUD = unique(LONGITUD), LATITUD = unique(LATITUD), .groups = "keep")

Dibujamos en el mapa de argentina los puntos que tenemos en el dataset agrupado. Para eso cargamos el dataset de departamentos que tiene info de toda la argentina

library(ggplot2)
departamentos <- st_read("data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp")

Reading layer `pxdptodatosok' from data source `/Users/antonellaschiavoni/Antonella/Maestria/EstadisticaEspacial/TPFinal/Datos/trabajo_final/smn/estadistica_espacial/data_departamentos/Codgeo_Pais_x_dpto_con_datos/pxdptodatosok.shp' using driver `ESRI Shapefile'
Simple feature collection with 527 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.02985 ymin: -90 xmax: -25.02314 ymax: -21.74506
Geodetic CRS:  WGS 84

#verificamos la proyeccion de departamentos
st_crs(departamentos)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

#como tiene otra proyeccion la pasamos a mercator
departamentos <- st_transform(departamentos, crs = 4326)

#ggplot() +
#  geom_sf(data = departamentos) +
#  geom_point(data = data_agg, aes(x = x, y = y), colour = "red", size = 1) +
#  theme_minimal() +
#  ggsave('plot.png',height = 25 , width = 30)

Convertimos data_agg a data.frame

df_data_agg = data.frame(data_agg)

Transformamos df_data_agg en un archivo geográfico utilizando el código de proyección mercator

data_sf = st_as_sf(df_data_agg, coords = c("LONGITUD", "LATITUD"), crs = 4326)

Validamos la clase del nuevo dataframe

class(data_sf)

[1] "sf"         "data.frame"

##TODO: graficar el mapa de densidades con los puntos de las estaciones en data.

GRAFICO INTERACTIVO

library(leaflet)

leaflet() %>%
  addTiles() %>%
  addCircles(data = data_sf, weight = 3)

NA

Analisis Estadistico de los datos

Analisis univariado

describe(data_agg$MEAN_VIENTO_KMH)

hist(data_agg$MEAN_VIENTO_KMH)

boxplot(data_agg$MEAN_VIENTO_KMH)

hist(data_agg$SD_VIENTO_KMH)

boxplot(data_agg$SD_VIENTO_KMH)

Analisis descriptivo de normalidad de los datos

hist(data_agg$MEAN_VIENTO_KMH)

boxplot(data_agg$MEAN_VIENTO_KMH)

qqnorm(data_agg$MEAN_VIENTO_KMH)
qqline(data_agg$MEAN_VIENTO_KMH, col=2)

hacemos el test de normalidad

shapiro.test(data_agg$MEAN_VIENTO_KMH)

    Shapiro-Wilk normality test

data:  data_agg$MEAN_VIENTO_KMH
W = 0.96947, p-value = 0.02215

El p-value obtenido es menor a 0.05, lo cual indica que los datos que tenemos no son normales

coordinates(data_agg) = ~LATITUD+LONGITUD
class(data_agg)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

plot(data_agg)

vt_exp = fit.variogram(v, vgm(125, "Exp", 30, 5))
vt_exp

plot(v , vt_exp)

vt_mat = fit.variogram(v, vgm(125, "Mat", 30, 5))
plot(v , vt_sph)

vt_exc = fit.variogram(v, vgm(125, "Exc", 30, 5))
plot(v , vt_sph)

attr(vt_exp, 'SSErr')

[1] 0.0126554

attr(vt_mat, 'SSErr')

[1] 0.0126554

attr(vt_exc, 'SSErr')

[1] 0.01653322